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We describe a method for analyzing the stochasticity in the non- stationary data for the beat-to-beat 
fluctuations in the heart rates of healthy subjects, as well as those with congestive heart failure. The 
method analyzes the returns time series of the data as a Markov process, and computes the Markov 
time scale, i.e., the time scale over which the data are a Markov process. We also construct an 
effective stochastic continuum equation for the return series. We show that the drift and diffusion 
coefficients, as well as the amplitude of the returns time series for healthy subjects are distinct from 
those with CHF. Thus, the method may potentially provide a diagnostic tool for distinguishing 
healthy subjects from those with congestive heart failure, as it can distinguish small differences 
between the data for the two classes of subjects in terms of well-defined and physically-motivated 
quantities. 
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Introduction 

Cardiac intcrbcat intervals fluctuate in a complex man- 
ner [1-7]. Recent studies reveal that under normal con- 
ditions, beat-to-beat fluctuations in the heart rate may 
display extended correlations of the type typically exhib- 
ited by dynamical systems far from equilibrium. It has 
been shown [2], for example, that the various stages of 
sleep may be characterized by long-range correlations in 
the heart rates, separated by a large number of beats. 

The analysis of the interbeat fluctuations in the heart 
rates belong to a much broader class of many natural, as 
well as man-made, phenomena that are characterized by 
a degree of stochasticity. Turbulent flows, fluctuations 
in the stock market prices, seismic recordings, the in- 
ternet traffic, pressure fluctuations in chemical reactors, 
and the surface roughness of many materials and rock [8] , 
are but a few examples of such phenomena and systems. 
A long standing problem has been the development of 
an effective reconstruction method for such phenomena. 
That is, given a set of data for certain characteristics 
of such phenomena (for example, the interbeat fluctua- 
tions in the heart rates), one would like to develop an 
effective equation that can reproduce the data with an 
accuracy comparable to the measured data. Although 
many methods have been suggested in the past, and con- 
siderable progress has been made, the problem remains, 
to a large extent, unsolved. 

In many cases the stochastic process to be analyzed 
is non- stationary. If the process also exhibits extended 
correlations, then deducing its statistical properties by 
the standard methods of analyzing such processes is very 
difficult. One approach to analyze such processes was 



proposed by Stanley and co-workers [1,3,5,20-24] and 
others [25-29]. They studied data for heart-rate fluctua- 
tions, for both healthy subjects and those with congestive 
heart failure (CHF), in terms of self-affine fractal distri- 
butions, such as the fractional Brownian motion (FBM). 
The FBM is a non-stationary stochastic process which in- 
duces long-range correlations, the successive increments 
of which are, however, stationary and follow a Gaussian 
distribution. The power spectrum of a FBM is given by, 
S(f) cx /~' 2ff+1 \ where H is the Hurst exponent that 
characterizes the type of the correlations that the data 
contain. Thus, one may distinguish healthy subjects from 
those with CHF in terms of the numerical value of H as- 
sociated with the data: negative or antipersistent corre- 
lations for H < 1/2, as opposed to positive or persistent 
correlations for H > 1/2. The analysis of Stanley and co- 
workers indicated that there may indeed be long-range 
correlations in heart-rate fluctuations data that can be 
characterized by the FBM and similar fractal distribu- 
tions. In addition, the data for healthy subjects seem to 
be characterized by H < 1/2, whereas those with CHF 
by H > 1/2. This was a significant discovery over the 
traditional methods of analyzing non-stationary data for 
heart-rate fluctuations. 

However, values of the Hurst exponent H associated 
with the two groups of subjects are non-universal. Thus, 
it would, for example, be difficult to distinguish the two 
groups of subjects if their associated Hurst exponents are 
both close to 1/2. In addition, the FBM is a non-self- 
averaging distribution, i.e., given a fixed Hurst exponent 
H, each realization of a FBM may be significantly dif- 
ferent from its other realizations with the same H. As 
a result, estimating H alone and characterizing the data 



by a FBM cannot enable one to predict the future trends 
of the data. One may also analyze such data by the de- 
tcrendcd fluctuating analysis [2-5] which, in many cases, 
is capable of yielding accurate and insightful information 
about the nature of the data. 

Recently, a novel method of analyzing stochastic pro- 
cesses was introduced [9-12]. It was shown that by an- 
alyzing stochastic phenomena as Markov processes and 
computing their Markov time (or length) scale (that is, 
the time scale over which the process can be thought of as 
Markov), one may reconstruct the original process with 
similar statistical properties by constructing an effective 
equation that governs the process. The constructed equa- 
tion helps one to understand the nature and properties 
of the stochastic process. The method utilizes a set of 
experimental data for a phenomenon which contains a 
degree of stochasticity, and constructs a simple equa- 
tion that governs the phenomenon [9-16]. The method 
is quite general; it is capable of providing a rational ex- 
planation for complex features of the phenomenon. More 
significantly, it requires no scaling feature. 

In this paper we describe a method for analyzing non- 
stationary data, and then utilize it to study the inter- 
beat fluctuations in the heart rates. We show that the 
application of the method to the analysis of interbeat 
fluctuations in the heart rates may potentially lead to 
a novel method for distinguishing healthy subjects from 
those with CHF. 

The plan of this paper is as follows. In the next section, 
we describe the method. We then utilize the method to 
analyze data for heart-rate fluctuations in human sub- 
jects. 

Markov Analysis of Non-Stationary Data 

Given a (discrete) non-stationary time series , we in- 
troduce a quantity x% , called the return of rj , defined by 

Xi = \n(r i+ i/ri) , (1) 

where is the value of the stochastic quantity at step i. 
If there are long-range positive correlations in the series, 
then Ti and r,+i are close in values and, therefore, we 
expect the series Xi to have very small values for all t. 
For white noise, as well as data that exhibit negative or 
anti-correlations, and ri+\ can be completely different 
and, therefore, the time series Xi will fluctuate strongly. 

Figures 1 and 2 present the typical data T{ and the 
corresponding returns Xi for healthy subjects and those 
with CHF. The number of data is of the order of 30,000- 
40,000, taken over a period of about 6 hours. It is evident 
that the returns series for the subjects with CHF has 
small amplitudes, implying that the T"j data set has long- 
range positive correlations, which is consistent with the 
previous analysis [1]. It can be verified straightforwardly 
that the series Xi is stationary, by measuring the stability 
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FIG. 1. Interbeats fluctuations of healthy subjects (top), 
and its returns (bottom). 
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FIG. 2. Interbeats fluctuations of subjects with congestive 
heart failure (top), and its returns (bottom). 

of its average and variance in a moving window (that 
is, over a period of time which varies over the length 
of the series). Due to the stationarity of the series 
x(t), we can construct an effective stochastic equation 
for the returns series of the two groups of subjects, and 
distinguish the data for healthy subjects from those with 
CHF. The procedure to do so involves two key steps: 

(1) Computing the Markov time scale (MTS) iju con- 
stitutes the first step. <m is the minimum time inter- 
val over which the data can be considered as a Markov 
process [9-12,17]. As is well-known, a given stochastic 
process with a degree of randomness may have a finite or 
even an infinite tjv/. To estimate the MTS 4m , we note 
that a complete characterization of the statistical proper- 
ties of stochastic fluctuations of a quantity x(t) in terms 
of a parameter t requires the evaluation of the joint prob- 
ability distribution function (PDF) P n (x\, t±; • • • ; x n ,t n ) 
for an arbitrary n, the number of the data points. If a 
stochastic phenomenon is a Markov process, an impor- 
tant simplification can be made as P n , the n-point joint 



PDF, is generated by the product of the conditional prob- 
abilities, tj+i \xi, ti), for i = 1, • • • ,n — 1. 

The simplest way to determine tu for stationary data 
is by using the least-square test. The rigorous mathe- 
matical definition of a Markov process is given [18] by 



P(xk,tk\xk-i,tk-V, • • • ; xi,ti;x , t ) 

= P(Xk,tk\Xk-l,tk-l) ■ 



(2) 



Intuitively, the physical interpretation of a Markov pro- 
cess is that it "forgets its past." In other words, only 
the closest "event" to Xk, say Xk-i at time tk-i, is rel- 
evant to the probability of the event Xk at tk- Hence, 
the ability for predicting the event Xk is not enhanced by 
knowing its values in steps prior to the the most recent 
one. Therefore, an important simplification that is made 
for a Markov process is that, the conditional multivariate 
joint PDF is written in terms of the products of simple 
two parameter conditional PDF's [18] as (3) 

P(xk,tk;xk-i,t k -i; ■ • ■;xi,t 1 \x ,t ) 



~[P{x i ,t i \x i _ 1 ,t i _ 1 ) 



(3) 



Testing Eq. (3) for large values of k is beyond the current 
computational capability. For k = 3 (three points or 
events), however, the working equation, given by, 



P(x 3 ,t 3 \x 2 ,t 2 ;xi,ti) = P(x 3 ,t 3 \x 2 ,t 2 ) , 



(4) 



should hold for any value of t 2 in the interval t\ < t 2 < t 3 . 
A process is then Markovian if Eq. (4) is satisfied for a 
certain time separation t 3 —t 2l in which case, tM = t 3 —t 2 . 
Thus, to compute the tM we use a fundamental theory of 
probability according to which we write any three-point 
PDF in terms of the conditional probability functions as, 

P(x 3 ,t 3 ;x 2 ,t 2 ;xi,t 1 ) 

= P{x 3l t 3 \x 2l t 2 \x ll t 1 )P{x 2l t 2 \x ll t 1 ). (5) 

Using the properties of Markov processes to substitute 
Eq. (5), we obtain, 

-P|VIarkov(a;3, ^3 ', X 2 , t 2 ; X\ , t \ ) 



= P{x 3 , t 3 \x 2 ,t 2 )P(x 2 ,t 2 \xi,ti). 



(G) 



We then compare the deviation of PMarkov from that 
given by Eq. (5). Using the least square method [10], 
we write: 

X 2 = / dx 3 dx 2 dx\X 



[P(x 3 ,t 3 ] X 2 , t 2 ] Xl,ti) - P]VIarkQv(x3, t 3 ] X 2 ,t 2 ] X\,ti)f 
2 i 2 
°" + ^Markov 

where a 2 and "Markov are ^ ne corresponding variances of 
terms in the nominator. Thus, one should plot the re- 
duced chi-square, xl = X 2 IN (with M being the number 
of degrees of freedom), as a function of the time scale 
t 3 — t 2 . Then, tm = t 3 — t 2 for that value of t 3 — t 2 for 
which xt either achieves a minimum or becomes flat and 
does not change anymore; see Figure 3. 

On the other hand, a necessary condition for a stochas- 
tic phenomenon to be a Markov process is that the 
Chapman-Kolmogorov (CK) equation (8), 

P(x 3l t 3 \xi,ti) = J dx 2 P(x 3l t 3 \x 2 ,t 2 ) P(x 2 ,t 2 \xi,ti) 

should hold for the time separation t 3 — t 2 , in which case, 
tM = t 3 — t 2 . Therefore, to test whether the time series 
x(t) is a Matkov process, one should check the valid- 
ity of the CK equation for describing the process using 
different x\ by comparing the directly-evaluated condi- 
tional probability distributions P (x 3 , t 3 \x-i, t\) with the 
one calculated according to right side of Eq. (8). 

(2) Estimation of the Kramers-Moyal coefficients is the 
second step of constructing an effective equation for de- 
scribing the series Xi. The CK equation is an evolution 
equation for the distribution function P(x,t) at any time 
t. When formulated in differential form, the CK equation 
yields the Kramers-Moyal (KM) expansion [18], given by, 



(7) 



(8) 
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p( X ,t) = y(~r[D^(x)p(x,t)] 



n=l 



O.T 



(9) 



The coefficients (x, t) are called the KM coefficients. 
They are estimated directly from the data, the condi- 
tional probability distributions, and the moments M^ n ' 
defined by, 

M<"' = -^ t J dx'(x' -x) n P(x',t + At\x,t), 

D(»)(a?,t) = 1 lim M {n) . (10) 
n\ At^o 

According to the Pawula's theorem, for a process with 
- all the with n > 3 vanish, in which case the 
KM expansion reduces to the Fokker-Planck equation, 
also known as the Kolomogrov equation [18]: 



P(x,t). (11) 



Here D^(x,t) is the drift coefficient, representing the 
deterministic part of the process, and D^> (x, t) is the 
diffusion coefficient that represents the stochastic part. 

We now apply the above method to the fluctuations 
in the human heartbeats of both healthy subjects and 




FIG. 3. xi values for a typical subject with CHF for several 
time scales. 

those with CHF. As mentioned in the Introduction, sev- 
eral studies [5,6,10-12,19-21] indicate that, under nor- 
mal conditions, the beat-to-beat fluctuations in the heart 
rate may display extended correlations of the type typ- 
ically exhibited by dynamical systems far from equilib- 
rium, and that the two groups of subjects may be distin- 
guished from one another by a Hurst exponent. We show 
that the drift and diffusion coefficients (as defined above) 
of the interbeat fluctuations of healthy subjects and pa- 
tients with CHF have distinct behavior, when analyzed 
by the method we propose in this paper, hence enabling 
one to distinguish the two groups of the subjects. 

We analyzed both daytime (12:00 pm to 18:00 pm) and 
nighttime (12:00 am to 6:00 am) heartbeat time series 
of healthy subjects, and the daytime records of patients 
with CHF. Our data base includes 10 healthy subjects 
(7 females and 3 males with ages between 20 and 50, 
and an average age of 34.3 years), and 12 subjects with 
CHF (3 females and 9 males with ages between 22 and 
71, and an average age of 60.8 years). Figures 1 and 2 
present the data. We first estimate the Markov time scale 
tM for the returns series of the interbeat fluctuations, 
using the chi-square method described above. In Figure 
3 the results for the xi values for a subject with CHF 
are shown. For the healthy subjects we find the average 
tjvf for the returns, for both the day- and nighttime data, 
to be (all the values are measured in units of the average 
time scale for the bcat-to-bcat times of each subject), 
tM = 10. On the other hand, for the daytime records 
of the patients with CHF, the estimated average t^j is, 
tM = 20. Therefore, the data for the healthy subjects are 
characterized by tM values that are smaller than that of 
the patients with CHF by a significant factor of 2. 

We then check the validity of the CK equation for sev- 
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FIG. 4. Test of Chapman-Kolmogorov equation for the 
time separation between ts and t\ equal to the Markov time 
scale, for x\ = —6 x 10 -2 , xi = 0, and x\ = 6 x 10~ 2 . Squares 
and triangles represent, respectively, the directly-evaluated 
PDF and that computed according to the right side of Eq. 
(8). For clarity, the PDFs are shifted in the vertical direc- 
tions. 

eral X\ triplets by comparing the directly-evaluated con- 
ditional probability distributions P(x3, ^Ixi, t\) with the 
ones calculated according to right side of Eq. (8). Here, 
x represents the returns. In Figure 4, the two differently- 
computed PDFs are compared. Assuming the statistical 
errors to be the square root of the number of events in 
each bin, we find that the two PDFs are statistically iden- 
tical. 

Using Eq. (10) directly we calculate the drift and dif- 
fusion coefficients, D^ 1 \x) and D^ 2 \x), for the entire set 
of data for the healthy subjects, as well as those with 
CHF. The corresponding D^'(x) and D^ 2 '{x) are dis- 
played in Figure 5. We find that, these coefficients pro- 
vide another important indicator for distinguishing the 
ill from the healthy subjects: The drift and the dif- 
fusion coefficients D^ 2 '{x) follow, respectively, linear and 
quadratic equations in x with distinct coefficients for the 
healthy subjects and patients with CHF. The analysis 
of the data yields the following estimate for the healthy 
subjects (averaged over the samples), 

D w {x) = -0.1a; , 

L> (2) (x) = 3.7 x 10~ 5 - 6.6 x 10~ 5 x + 0.06a; 2 , (12) 

with —0.15 < x < 0.15, whereas for the patients with 
CHF we find that, 

D (1) (x) = -0.06x , 

L> (2) (a:) = 8.6 x 10~ 6 - 2.7 x 10~ 5 x + 0.03x 2 . (13) 
with -0.04 < x < 0.04. 



We find two important differences between the heart- 
beat dynamics of the two classes of subjects: 



(1) Compared with the healthy subjects, the drift and 
diffusion coefficients for the patients with CHF are small. 

(2) The fluctuations of the returns for healthy subjects 
are distinct from those with CHF. They also fluctuate 
over different intervals, indicating that the returns data 
for the healthy subjects fluctuate over large interval. The 
fluctuations intervals are, —0.04 < x < 0.04 and —0.15 < 
x < 0.15 for patients with CHF and healthy subjects, 
respectively. Hence, we suggest that one may use the 
drift and diffusion coefficients magnitudes, as well as the 
fluctuations intervals for the returns, for characterizing 
the dynamics of human heartbeats, and to distinguish 
healthy subjects from those with CHF. 
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I. DISCUSSIONS 



Lin [30] argued that the daytime heart rate variability 
of healthy subjects may exhibit discrete scale-invariancc 
(DSI). A stochastic process x(t) possesses continuous 
scale-invariant symmetry if its distribution is preserved 
under a change of variables, t — > Xt and x — > x/fi, where 
A and \x are real numbers, so that, 

x{t) = -x(Xt) . (14) 
A* 

If Eq.(14) holds only for a countable (discrete) set of val- 
ues of A, x(t) is said to possess DSI, which implies a 
power-law behavior for x(t) that has a log-periodic cor- 
rection of frequency 1/logA, so that 

x(t)=tiF(logt/log\), (15) 

with, 7 = log /i/ log A, with F(x) = F(x + 1) being a pe- 
riod scaling function. Generally speaking, one may write, 
x(t) = c(t)t^, with, £ = 7 + 2n7ri/logA, with n = 1,2, ■ • ■ 
The existence of log-periodicity was first suggested by 
Novikov [31] in small-scale energy cascade of turbulent 
flows. It has been argued [32] that log-periodicity may 
exist in the dynamics of stock market crashes [33] , turbu- 
lence [34] , earthquakes [35] , diffusion in disordered mate- 
rials [36,37], and in fracture of materials near the macro- 
scopic fracture point [38] . The log-periodicity, if it exists 
in the heart rate variability (HRV) , implies the existence 
of a cascade for the multifractal spectrum of HRV, pre- 
viously reported by others. However, Lin's method, nei- 
ther provides a technique for distinguishing the HRV of 
healthy people from those with CHF, nor can it predict 
the future behavior of HRV based on some data at earlier 
times. 

The method proposed in the present paper is different 
from such analyses in that, the returns for the data are 
analyzed in terms of Markov processes. Our analysis does 
indicate the existence of correlations in the return which 
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FIG. 5. The drift and diffusion coefficients, D^'(x) and 
D^ 2 \x), estimated by Eq. (8). For the healthy subjects (tri- 
angles) and for patients with CHF (squares), D^(x) and 
D' 2 ' (x) follow linear and quadratic equations in x. 

can be quite extended (and is characterized by the value 
of the Markov time scale £m)- 

II. SUMMARY 

We distinguish the healthy subjects from those with 
CHF in terms of the differences between the drift and 
diffusion coefficients of the Fokker-Plank equations that 
we construct for the returns data which, in our view, pro- 
vide a clearer and more physical way of understanding 
the differences between the two groups of the subjects. 
In addition, the reconstruction method suggested in this 
paper enables one to predict the future trends in the re- 
turns (and, hence, in the original series r^) over time 
scales that are of the order of the Markov time scale tu- 



None of the previous approaches for analyzing the data 
could provide such a reconstruction method. 

We also believe that, the computational method that 
is described in this paper is more sensitive to small dif- 
ferences between the data for healthy subjects and those 
with CHF. As such, it might eventually provide a diag- 
nostic tool for detection of CHF in patients with small 
amounts of data and in its initial stages of development. 
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